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Abstract 

We compute the low-lying spectrum of the staggered Dirac operator above and below the finite 
temperature phase transition in both quenched QCD and in dynamical four flavor QCD. In both 
cases we find, in the high temperature phase, a density with close to square root behavior, p(A) ~ 
(A — Ao) 1 / 2 . In the quenched simulations we find, in addition, a volume independent tail of small 
eigenvalues extending down to zero. In the dynamical simulations we also find a tail, decreasing 
with decreasing mass, at the small end of the spectrum. However, the tail falls off quite quickly and 
does not seem to extend to zero at these couplings. We find that the distribution of the smallest 
Dirac operator eigenvalues provides an efficient observable for an accurate determination of the 
location of the chiral phase transition, as first suggested by Jackson and Verbaarschot. 
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1 Introduction 



The correlations of Dirac operator eigenvalues in QCD and related theories have been shown to have a 
fascinating relation to Random Matrix Theory. There are two very different domains of interest here. 
One is the so-called "bulk" of the eigenvalue spectrum of the Dirac operator, far from both the infrared 
and the ultraviolet ends. The other is the so-called "hard edge" at A ~ 0, i.e. the infrared end of the 
spectrum in theories with spontaneous breaking of chiral symmetry. The relevance of Random Matrix 
Theory in describing eigenvalue correlations of the Dirac operator in the bulk of the spectrum was first 
demonstrated by Halasz and Verbaarschot Q, and it has since been confirmed by numerous lattice 
gauge theory studies 0. From a theoretical point of view, these results in the bulk of the spectrum 
remain to be better understood. In the other case, near A ~ 0, the connection between universal 
Random Matrix Theory results and the QCD Dirac operator spectrum ||, |], |5| is by now firmly 
established, and has also been extensively checked in lattice gauge theory simulations ||. Already 
in the original work of Shuryak and Verbaarschot j|] it was shown that the pertinent chiral Random 
Matrix Theory partition function coincides exactly with the effective field theory partition function in 
the relevant microscopic limit. More recently an explicit relationship between the universal Random 
Matrix Theory eigenvalue distributions and those of the QCD Dirac operator has been established 
|7|, ||, the precise link being the partially quenched chiral condensate ||, [|. In this domain Random 
Matrix Theory is an intriguing alternative description of exactly the same phenomena that can be 
derived from the effective QCD partition function in the large- volume limit where V <C |jTC|| . 

There exists also an interesting physical situation that forces us to reconsider the two different do- 
mains of the Dirac operator spectrum simultaneously. This is near the finite-temperature chiral phase 
transition, where the analytical connection between the effective QCD partition function and Ran- 
dom Matrix Theory breaks down even for the part of the spectrum that is close to A = 0. This 
situation is readily confronted in lattice gauge theory. Indeed, it is for staggered fermions even rig- 



orously proven [11], that chiral symmetry is restored at high temperature. The low- lying spectrum 
of the Dirac operator must therefore be quite different in the high temperature phase as compared 
with zero temperature. In particular, the absence of chiral symmetry breaking implies, through the 
Banks-Casher relation, 

{W) = 7T P (0) , (1) 

that the density of eigenvalues at zero vanishes. This could happen either just at that point of A = 
(as, for example, in the free theory where p(X) ~ A 3 ), or the spectrum could develop a gap. The chiral 
phase transition occurs at the temperature T c where the density of Dirac operator eigenvalues just 
reaches zero at A = 0. If the transition is continuous, this will happen smoothly as the temperature 
T is increased towards T c . In such a case, an important question to settle is the precise power-law 
behavior of the spectral density of the Dirac operator right at T = T c , because this can be related to 
the critical exponents of the phase transition [^] . The same chiral Random Matrix Theory that yields 
universal microscopic spectral correlators which exactly coincide with those of the Dirac operator at 
T = can be tuned in such a way as to just reach (multi-critical) points where [13] 



p(X) ~ X 2k (2) 

near A = 0. Here k is an integer labeling the multi-criticality. At such points universality of microscopic 
spectral correlators still holds in the Random Matrix Theory context, but there is no justification for 
assuming that these results are relevant for the Dirac operator spectrum at T = T c .[] There is also a 

1 In particular, this behavior is only compatible with a continuous phase transition. But at a more fundamental 
level, there is simply no longer any relation between the chiral Random Matrix Theory ensemble and the effective QCD 
partition function for temperatures T that do not satisfy T <C Aqcd ■ 
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schematic Random Matrix model for the finite-temperature behavior of the Dirac operator spectrum 
[14]. It gives a different behavior at T = T c : p(X) ~ A 1 / 3 . Also here there is no physical justification for 



using it in connection with the Dirac operator spectrum at finite temperature, but it is an interesting 
model that depends on just one deterministic external parameter, and we shall therefore return to it 
below. 

Suppose, for a moment, that the Dirac spectrum actually develops a gap around A = above T c . 
In Random Matrix Theory the end of a spectrum around such a gap is referred to as a "soft edge" . 
Generally, the (macroscopic) density of eigenvalues near a soft edge behaves as (for A > Aq) p5fl: 



p(X) oc (A-A ) 2m+1 / 2 (3) 

with Ao being the location of the edge. Here m is an integer that labels universality classes classified 
by their Random Matrix Theory potentials (m parameters in the potentials must be tuned in order to 
reach each class). Thus the generic behavior, without any fine tuning, corresponds to m = 0, which 
gives a square root approach to the soft edge: 

p(A) oc (A-A ) 1/2 (4) 

Random Matrix Theory actually gives a more detailed, microscopic, description. This arises from a 
blowing-up of the eigenvalue density function around the soft edge Ao with a rescaling according to the 
macroscopic behavior (|J). For example, for the generic m = universality class, the corresponding 
microscopic eigenvalue density is 

p(X) oc X[Ai(-X)} 2 + [Ai'(-X)] 2 , (5) 

where X oc (A — Ao)A^ 2 / 3 and Ai(x) is the standard Airy function. Here N denotes the size of 
the matrix, and the rescaling by N 2 ^ 3 is required in order to spread out the increasing number of 
eigenvalues to obtain one well-defined limiting function. From the known asymptotic behavior of the 
Airy functions one finds: 

^ + O (x) for X ^ oo 



o(X) ~ < * \ x > (6) 

\_^|X|- 1 /2 exp (_4|x[ 3 / 2 /3) forX^-oo 

Thus, at the microscopic level the spectrum is not cut off sharply at Ao, but has an exponentially 
suppressed tail beyond. Further, the square root behavior of the eigenvalue density is modulated by 
wiggles corresponding to the distribution of particular eigenvalues (smallest, second, third, etc.). (See, 
for example, Fig. ||.) 

Until a few months ago, there existed only one low-statistics investigation of the low-lying Dirac 
eigenvalue spectrum for staggered fermions in the high temperature phase [17]. It did not unequivocally 



establish the existence of a gap. For example, some low eigenvalues were found. However, they could 
possibly be attributed to "would-be zero modes" from global gauge field topology, shifted away from 



zero by the explicit chiral symmetry breaking of staggered fermions at finite lattice spacing [18]. This 
is a general problem with staggered fermions that also we must face here: at finite lattice spacing, 
the index theorem is not valid for staggered fermions, and gauge field topology (whichever way one 
defines it on a discrete lattice) does not give rise to exact zero modes of the staggered Dirac operator. 
At T = no trace of non-trivial gauge field topology on the lowest-lying spectrum of staggered 



eigenvalues has been found [19], except in the 2-d Schwinger model at fairly small lattice spacing |2Cj. 
Finite temperature, which causes a depletion of genuine non-zero eigenvalues near A ~ 0, further 
complicates this issue of a mix-up with would-be zero modes. 
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Very recently a study of just the low-lying Dirac operator eigenvalues near T c has actually indicated 



the presence of a gap at large lattice volumes [21 1. Again the statistics was rather limited, and only 



a small number of the lowest-lying eigenvalues could be included (varying between 8 and 10). The 
present study will in many ways follow the same lines of thought as Ref . |2l]] , but we shall have much 
larger statistics, and we shall also probe some different aspects. In particular, we are also interested 
in seeing whether quenching causes a different behavior of the smallest eigenvalues near T c compared 
with dynamical fermions. 

Last year a study of low-lying eigenvalues of the overlap Dirac operator in quenched finite temperature 
gauge theories appeared |22| |. Overlap fermions are well suited for such an investigation since they do 
not suffer from explicit chiral symmetry breaking even at finite lattice spacing and since they have exact 
zero modes in topologically non-trivial gauge fields. Ref. J22[ | found that effects of topology persisted 
in the high temperature phase, although strongly suppressed compared to the low temperature phase. 
More interestingly, an accumulation of low eigenvalues with an apparently finite density in the infinite- 
volume limit was found very near A ~ even above the (quenched) phase transition temperature T c . 
The statistical properties of the smallest group of eigenvalues were consistent with them being due 



to a dilute gas of instantons and anti-instantons [22|. These results have led to the speculation that 



chiral symmetry might remain broken even in the high temperature phase of quenched QCD with 
overlap fermions. 

In this paper we describe high statistics investigations of the low-lying eigenvalue spectrum of the 
staggered Dirac operator for SU(3) gauge group at finite temperature. We do this for both quenched 
(section ||) and dynamical QCD with four flavors of staggered fermions (section ||). The interest in 
the quenched case lies primarily in checking whether also staggered fermions, although insensitive to 
topology at the gauge couplings we can investigate here, show signs of unusual behavior of the smallest 
Dirac eigenvalues above T c (as was the case for overlap fermions [22|). Both quenched and unquenched 



simulations give rise to Dirac operator spectra that can be compared with Random Matrix Theory. In 
particular, if the Dirac operator spectrum exhibits a gap, is it of the soft-edge kind of Random Matrix 
Theory? Are there indications that Random Matrix Theory provides a more accurate description of 
the Dirac operator spectrum as the three- volume is increased? We shall try to answer these questions 
in what follows. 



2 Quenched QCD at high temperature 

For our quenched finite temperature Monte Carlo simulations we have used lattices with temporal 
extent Nt = 4 and up to three spatial volumes: 8 3 , 12 3 and 16 3 . For JVj = 4 the deconfinement phase 
transition for SU(3) pure gauge theory with Wilson action has been very accurately determined, 
occurring at lattice gauge coupling C = 5.6925(2) p3| . It has for long been assumed that the chiral 
phase transition of the quenched theory occurs at exactly this deconfinement phase transition point 
of the pure gauge theory, but this has recently been challenged p^|. 

We now give some technical details of our simulations. The gauge field configurations were generated 
with a mixture of overrelaxation and heat bath updates, and were analyzed after every 20-th heat 
bath sweep. On each configuration we computed the 50 lowest-lying eigenvalues, and in some cases 



even more, using the variational Ritz functional method [24|. Many ensembles consisted of several 
thousand configurations. Gauge coupling, lattice size and statistics of our ensembles is summarized 
m Table g. 
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V 


P 


^configurations 


8 3 x 4 


5.5 


5931 


8 3 x 4 


5.66 


5544 


8 3 x 4 


5.695 


1769 


8 3 x 4 


5.71 


1145 


8 3 x 4 


5.72 


1426 


8 3 x 4 


5.73 


2021 


8 3 x 4 


5.75 


7029 


12 3 x 4 


5.75 


1933 


16 3 x 4 


5.75 


692 


8 3 x 4 


5.8 


5335 


8 3 x 4 


5.85 


6149 


8 3 x 4 


5.9 


4263 


12 3 x 4 


5.9 


1755 


16 3 x 4 


5.9 


824 



Table 1: Details of our quenched ensembles. 
As is well-known, the staggered Dirac operator 

= $e,o+l)o,e (7) 

is anti-hermitian, with purely imaginary eigenvalues iX that come in pairs of opposite sign. In eq. (0) 

ty(x) = {-l)^<» Xv (8) 
are the usual phase factors for staggered fermions. Denoting 

c(x) = (-1)£,*» , (9) 

we have also explicitly shown howl) connects even sites, i.e. those with e{x) = +1, with odd sites, those 
with e(x) = —1, and vice versa. This means that the operator — L/) 2 is hermitian and positive semi- 
definite. The sign function e(x) defined above plays the role of 75 in the continuum: it anticommutes 
withi/f : {I), e} = 0. As — L/> 2 does not mix between even and odd lattice sites, we need only compute 
the eigenvalues, on, say, the even sublattice. One easily sees that if ip e is a normalized eigenvector of 
—I]) 1 with eigenvalue A 2 , then ip = j0 Ote ij} e is a normalized eigenvector of —L/> 2 with eigenvalue A 2 , 
and non-zero only on odd sites. Moreover, as we will never encounter exact zero modes on genuine 
quantum configurations, there is no difficulty with the above definition of ip Q . We of course make use 
of these properties, and hence compute only the (positive) eigenvalues of —I) 2 restricted to the even 
sublattice, and then take the (positive) square root. All eigenvalues to be shown in the following thus 
have an equal number of negative companions, of the exact same magnitude. 

The spectral density of the Dirac operator is given by 

pW = ^<E 5 ( A - A ")> . ( 10 ) 
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Figure 1: The measured distribution of the 50 lowest eigenvalues in quenched QCD at finite tempera- 
ture for several /3's on 8 3 x 4 lattices. One sees a clear suppression of the spectral density at the origin 
jo(0) as the temperature (coupling) is increased, indicating a diminishing chiral condensate according 
to the Banks-Casher formula. For large temperatures the spectrum is heavily suppressed near the 
origin, although in this quenched case a small tail always seems to extend towards A = (see below). 
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Figure 2: The density of the smallest eigenvalues in the high temperature phase of quenched QCD, at 
(3 = 5.9, for spatial volumes 8 3 , 12 3 and 16 3 . In the two large volumes we measured 100 eigenvalues 
and only 50 in the small volume. Since the eigenvalue density increases with increasing volume, the 
distributions are correspondingly cut off at smaller A. There is no observable change in the eigenvalue 
density as the three- volume is increased. This holds even at the very small tail extending towards 
A = 0, as shown in the magnified plot inserted. 

and it is readily computed numerically from our Monte Carlo simulations by a binning of the measured 
eigenvalues per configuration at convenient small intervals. 

In Fig. [l| we show the density of low lying eigenvalues on 8 3 x 4 lattices for several j3 values between 
5.5 (in the confined phase) to 5.9 (in the deconfined phase). In each case we computed the 50 lowest 
positive eigenvalues i\ of the staggered Dirac operator. The plots are normalized by the following 
condition, / °° p(X)d\ = ^eigenvalues/^. In the confined phase it is quite evident that the density at 
zero would be non-zero in the thermodynamic limit. As the temperature is increased, the density at 
zero decreases. There is a qualitative change of the eigenvalue density once the temperature is increased 
above T c . Beyond the first few eigenvalues the eigenvalue density assumes a shape compatible with 
a square root behavior (||). However, a sizable tail, decreasing with increasing temperature, of small 
eigenvalues persists, which seems to extend all the way down to zero. 

In Fig. H we compare the eigenvalue density in the high temperature phase, at (3 = 5.9, for several 
spatial volumes. As can be seen the eigenvalue density is volume independent to a surprising degree. 
In particular, the tail of small eigenvalues is seen to be volume independent and appears to persist in 
the thermodynamic limit. We note that the tail is much larger than would be compatible with the 
exponential tail from the Airy function behavior of random matrix theory at a soft edge eq. (|6|) (see 
Fig. ^ in section 3 below for an example). 

It is tempting to speculate that the tail seen here in the quenched case with staggered fermions is a 
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Figure 3: Distributions of the unfolded level spacings between the first few eigenvalues for (3 = 5.75 
on an 8 3 x 4 lattice, compared to the Wigner surmise. The agreement with Random Matrix Theory 
is clear once we go beyond the first few eigenvalues. However, there is not very good agreement for 
the first 2-3 eigenvalues at this volume. 



reflection of the accumulation of small eigenvalues seen with overlap fermions in |22[ and attributed 
there to a dilute gas of instantons and anti-instantons. In the staggered case, the explicit chiral 
symmetry breaking at finite lattice spacing shifts the small modes, presumably resulting in the observed 
tail. 



We also looked at the unfolded level spacing between individual eigenvalues i and i + 1 

s _ Aj + i — Xj 
(Aj+i — Aj) 

In Fig. H| the distribution of the Sj between eigenvalue i and i + 1, i = 1..6 is shown and compared to 
the expected distribution, the Wigner surmise 

P( s ) = . (12) 

Evidently, the level spacings between the first eigenvalues are not very accurately given by Random 
Matrix Theory correlations. But as we move into the bulk, say for i > 5 in the case of the 8 3 x 4 
lattice, the agreement with Random Matrix Theory becomes almost perfect. It should be stressed 
here that this is a volume-dependent statement. For larger volumes one needs to go beyond more 
eigenvalues starting at the soft edge before one finds good agreement. This is consistent with the fact 
that there is a definite scale around the soft edge, where eigenvalue correlations are poorly described 
by Random Matrix Theory. Going to larger volumes simply forces more eigenvalues into that region. 
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Figure 4: The lowest eigenvalue, averaged over 5 configurations, on an 8 3 x 4 lattice with am q = 0.025 
as a function of (3. One sees a clear change in behavior around (3 ~ 5.02, indicating the restoration of 
chiral symmetry around this temperature. Just this first Dirac operator eigenvalue can thus serve as 
an excellent indicator of the chiral phase transition point. 

3 QCD with four dynamical fermions at high temperature 

For our dynamical simulations we chose to work with n/ = 4 flavors. Four is the "natural" number 
of flavors for staggered fermions in the continuum limit, and an efficient exact simulation algorithm 
can be used, the Hybrid Monte Carlo algorithm. In addition, the four flavor theory is known to have 
a rather strong first order finite temperature phase transition p5] ]. Therefore, tunnelings into the low 
temperature phase are strongly suppressed already on rather small systems and at temperatures quite 
close to T c . 

We made dynamical simulations with quark masses am q = 0.1, 0.05, 0.025, 0.01 and 0.002, and for 
couplings (3 in the high temperature phase (here o is the lattice spacing). Interestingly, the lowest 
eigenvalue provides a sensitive method for determining the critical coupling f3 c (am q ): in Fig || the 
lowest eigenvalue on a 8 3 x 4 lattice with am q = 0.025 is plotted against (3. (3 ranges from 4.96 to 
5.07 in intervals of 0.002. Each point is an average over 5 configurations at that /3-value. A rise is 
observed between [3 = 5.012 and (3 = 5.022, which can be clearly interpreted as the chiral symmetry 
restoring finite temperature phase transition. (For published values of the critical couplings (3 c {am q )) 
see Ref. f25|.) The possibility of using the magnitude of the smallest Dirac operator eigenvalues as 
probes for chiral symmetry restoration was suggested by Jackson and Verbaarschot [14| on the basis 



of a similarly observed behavior in a Random Matrix Theory context that we will also discuss below. 

Most of our analysis with dynamical fermions was performed at (3 = 5.2, which is above (3 C for all 
values of m q we used. Typically, we analyzed ~ 3000 configurations for each volume and m q , extracting 
the 50 smallest eigenvalues. Ensemble details are summarized in Table [2|. Note the large statistics 



S 



V 


P 


am q 


^configurations 


8 3 x 4 


5.1 


0.1 


8897 


8 3 x 4 


5.1 


0.025 


2909 


8 3 x 4 


5.1 


0.01 


4050 


8 3 x 4 


5.2 


0.1 


4348 


8 3 x 4 


5.2 


0.05 


2149 


8 3 x 4 


5.2 


0.025 


29632 


12 3 x 4 


5.2 


0.025 


13929 


8 3 x 4 


5.2 


0.01 


4050 


12 3 x 4 


5.2 


0.01 


2882 


8 3 x 4 


5.2 


0.008 


4200 


8 3 x 4 


5.2 


0.005 


7948 


8 3 x 4 


5.2 


0.002 


6950 


8 3 x 4 


5.3 


0.025 


4549 


8 3 x 4 


5.4 


0.025 


3746 



Table 2: Details of our ensembles with dynamical fermions. 
gathered for (3 = 5.2 and am q = 0.025. 

As in the quenched case we find a small tail at the lower end of the eigenvalue distribution, reaching 
towards A = from the main bulk of the distribution. This tail also appears to be volume independent 
(see Fig. ^). However, in this case it is somewhat more suppressed than in the quenched case, and for 
small values of the quark mass it does not extend all the way to A = 0, see Fig. ||. 

In Fig. [7| we show eigenvalue distributions measured at j3 = 5.2 and various am q . Note that for 
larger am q , j3 c is also larger, which causes a shift of the whole spectrum as am q is varied. Therefore, a 
direct quantitative comparison of the distributions in Fig. fj] is not straightforward. However, we note 
that the distributions at am q = 0.025 and am q = 0.002 are almost on top of each other, indicating 
that these distributions are very close to the am q = distribution. 

Furthermore, as in the quenched case described in the previous section, we find the "macroscopic" 
behavior of the eigenvalue density seemingly compatible with a square root form. By fitting eq. (||) to 
the bulk of the spectrum (leaving out the tail), we obtain different powers, depending on how much 
of the tail we choose to cut off. This fitting has been done for V = 8 3 x 4 at (3 = 5.2 and am q = 0.025 
and the resulting values can be seen in Table ||. There is a clear tendency towards a slightly higher 
power (= 2m + 1/2) than the square root which would be expected from Random Matrix Theory. In 
Fig. |H we show the fit with the cut at 0.1. For comparison, we include the distribution of the first 
eigenvalue in the figure. Clearly, the tail at small A is caused by the excessive width of the distribution 
of the smallest eigenvalue. 

In Fig. [9| we compare the spectral density for quark mass am q = 0.025 with the prediction (||) of the 
Random Matrix Theory for the density near a soft edge. In order to make the comparison possible, 
we rescale the spectral density as 

where we determine Aq and K from a fit to a square root: J ^-Ky/ (A — Aq). Here V is the lattice vol- 
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Figure 5: Quenched spectrum of f3 = 5.9 and a dynamical spectrum at f3 = 5.4 with am q = 0.025, 
both on an 8 3 x 4 lattice. A blowup of the y-axis near the end of the tail is also shown. There is 
a much stronger suppression of small eigenvalues in the case of dynamical fermions, perhaps due to 
fewer would-be zero modes in topological non-trivial gauge field configurations. 



Cut at 
0.09 

0.095 
0.1 

0.105 
0.11 

0.115 
0.12 

0.125 
0.13 
0.14 



m 

0.055 ±0.003 
0.049 ±0.004 
0.040 ±0.005 
0.036 ±0.006 
0.041 ±0.008 
0.035 ±0.010 
0.043 ±0.013 
0.040 ±0.016 
0.033 ±0.021 
0.027 ±0.034 



2m + 1/2 
0.609 ±0.006 
0.598 ±0.007 
0.580 ±0.010 
0.572 ±0.011 
0.582 ±0.016 
0.570 ±0.019 
0.585 ±0.027 
0.581 ±0.032 
0.567 ±0.041 
0.554 ±0.068 



Table 3: The dependence of the fitted power on the A where the small-A tail is cut off. The data are 
from 8 3 x 4, (3 = 5.2 and am q = 0.025 lattices, using the 50 smallest eigenvalues. The fits are almost 
compatible with a square-root behavior of the spectrum at this soft edge, but there is a consistent 
small upward shift in the exponent. 
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Figure 6: The density of the 50 smallest eigenvalues in the high temperature phase, (3 = 5.2 and 
arriq = 0.025 for spatial volumes 8 3 and 12 3 . The normalization of the distributions is as in Fig. ||. As 
in the quenched simulations, we observe no significant change in the eigenvalue spectrum at all as the 
three- volume is increased. 
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Figure 7: The density of the 50 smallest eigenvalues at (5 = 5.2 for quark masses am q = 
0.1,0.05,0.025,0.002 on 8 3 x 4 lattices. While there is some shift with decreasing quark mass for 
the larger values of am q , the spectral density seems to approach a limiting function, which we can 
identify as the chiral limit of the spectrum. Note that the large-mass behavior of the spectrum is al- 
most identical to the quenched spectrum in this region, in agreement with the notion that on this scale 
of eigenvalues the large-mass fermion has completely decoupled, and effectively become quenched. 
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Figure 8: Power fit to the spectral density for = 5.2 with am q = 0.025 on an 8 3 x 4 lattice with cut 
at 0.1 (see Table |3|). Also included is the distribution of the first eigenvalue. The distribution of this 
smallest eigenvalue is very different from that demanded by the Airy- function behavior (see below). 



ume. We see that the tail predicted by the Random Matrix Theory is much more strongly suppressed 
than the measured distribution. Moreover, we can also observe that the density of the eigenvalues 
themselves does not match: the measured distribution includes 50 eigenvalues, whereas the function 
(||) has only 40 'wiggles' in the A-range of the distribution. However, if we attempt to perform the 
comparison by matching the number of eigenvalues /wiggles, the overall fit becomes worse. In itself, 
this mismatch is not surprising, when we remember that the overall shape is not well described by a 
square root behavior in the first place (see Table |3|) . 

In the quenched simulations we saw that correlations between the first 5 or 6 eigenvalues were not 
very accurately given by Random Matrix Theory. In Fig. ^ the distribution of Si, defined in eq. (11), 
is shown for a dynamical simulation with (5 = 5.2 , am q = 0.025 on a 8 3 x 4 lattice. Here we see, 
on the same lattice volume, a clear deviation from the Wigner surmise only for the first two or three 
spacings, s%, S2 and S3. Again, comparisons can only be made at equal volumes, as larger volumes 
imply more eigenvalues in the region around the soft edge where correlations are poorly described by 
Random Matrix Theory. 

There is an interesting physical consequence of a genuine gap in the Dirac eigenvalue spectrum. 
Consider the difference between the (it) susceptibility \P an d the susceptibility of its scalar partner 
(ao), which we denote by xs WA - I n a manner similar to the Banks-Casher formula for the chiral 
condensate one finds that this difference can be written 
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Figure 9: Comparison of the spectral density for (5 = 5.2, am q = 0.025 on a 8 3 x 4 lattice, with the 
Random Matrix Theory prediction for a soft edge, as in eq. (pi). 
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Figure 10: Distributions of the unfolded level spacings between the first few eigenvalues in the dynam- 
ical case, compared to the Wigner surmise. As in the quenched case, we find a small disagreement for 
the first eigenvalues. But in contrast to the quenched case the agreement sets in after almost just the 
2nd eigenvalue, at the same volume. 
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where the spectral density p(A; m) includes the zero modes due to topology as well. In the infinite- 
volume limit this contribution from exact zero modes should vanish, and we should be left with the 
integral over non-zero modes. If the spectral density has a gap around the origin this means that u 
vanishes in the chiral limit. Because axial U(l) rotates xp Xs an d vice versa, this would imply 
a restoration of axial U(l) at these high temperatures [26, 27]. Conversely, if one believes that this 
axial U(l) symmetry is not restored at high temperature, one gets constraints on the behavior of 
the spectral density of the Dirac operator near the origin. For a power-law behavior near A ~ a 
non-vanishing u> in the infinite-volume chiral limit can only be supported for p(A, 0) ~ A" with a < 1, 
and in fact uj would only remain finite in that limit if a = 1. 



Of course, the fact that we appear to see a gap in the eigenvalue spectrum with staggered fermions at 
this particular coupling does not imply that this gap persists in physical units as we take the continuum 
limit. If not, we are here studying a pure lattice artifact. The only way to test this is to study the shift 
in the apparent cut-off eigenvalue Ao as we change the lattice spacing (this is, however, far beyond 
what we can do at present). A very different uncertainty comes from the fact that we are working 
with fermions that are almost insensitive to topology at the couplings and volumes available to us. 
This means that there are some would-be zero modes mixed up with our regular eigenvalues near the 
origin. At zero temperature we found that these would-be zero modes somewhat surprisingly behaved 
as the non-zero eigenvalues at realistic couplings and volumes j0|] . But there is no guarantee that this 
is the case at finite temperature. This means that the small tail extending to zero in our quenched 
simulations, and some of the smallest eigenvalues in the simulations with dynamical fermions may be 
due to these would-be zero modes. In particular, the disagreement with Airy-function behavior very 
close to the soft edge may be partly due to these would-be zero modes. 



4 Random Matrix Theory 

We have just shown our lattice gauge theory data for the spectrum of the Dirac operator below, around, 
and above T c , and compared it with some of the analytical formulas of Random Matrix Theory in the 
limit where the size of these matrices N goes to infinity. It is of interest to see how a simple model of 
a chiral phase transition in Random Matrix Theory behaves at finite N, as some external parameter 
(mimicking temperature T) is changed. The model we shall focus on here was proposed by Jackson 



and Verbaarschot in the first of Ref. 14] and has also been studied in Ref. [28]. It is based on a 



modified chiral ensemble of N x N complex matrices W, with partition function 

N f 



in a sector of topological charge zero. Here 



JdW J| det(M - im/) exp \-NTr(WW ] )\ (15) 



W* + T 



M = [W + T ) ™ 

is a 2A^ x 2N block hermitian matrix, and the external deterministic parameter T is playing a role 
reminiscent of temperature. In the normalization chosen here, a continuous "phase transition" occurs 
at T c = 1, where the spectral density p{\) of the eigenvalues of the random matrices just reaches zero 
flU] . For T > T c the spectrum becomes two-banded, with a gap surrounding zero. The global shape of 
p(A) in this model will be very different from the actual (macroscopic) Dirac operator spectral density. 
But the interesting feature lies in having here a simple model which qualitatively seems to describe 
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Figure 11: Full spectra of 200 x 200 matrices at T 
gap indeed develops around A = at T = 1. 



0, 0.5, 0.75, 1.0, 1.2, 2.0. In this simple model a 



some of the observed behavior of the staggered Dirac operator with Nf = 4, in particular the apparent 
gap in the density for T > T C .Q 

It is very simple to do quenched (i.e, Nf = 0) numerical simulations of the above Random Matrix 
model, as it just corresponds to generating an ensemble of complex matrices with Gaussian weight. 
Such matrices are "maximally random" in that their matrix elements are independently of Gaussian 
distribution. This allows us to choose very large random matrices numerically, and then finding the 
eigenvalues of the hermitian matrix M. In this way we have studied the detailed behavior of the 
finite- iV Random Matrix model (15) just below T c , at T c , and just above T c (where the two bands 
separate as the gap develops) .0 



We show some of our numerical results in Figs. 11 and 12. First, in Fig. 11 we display a se- 
quence of the macroscopic spectral density with increasing magnitude of the parameter T: T = 
0,0.5,0.75,1.0,1.2,2.0. The plots were made by diagonalizing 10000 200 x 200 matrices. At T = 
this macroscopic Random Matrix Theory spectrum is just the Wigner semi-circle law (we display only 
the A > part): 

p(A,T = 0) 



2vr 



A 2 



(17) 



As T increases, a dip in the spectral density around A = slowly develops, and it subsequently 
turns into a gap. Of course, this macroscopic Random Matrix Theory spectrum is totally unlike the 
macroscopic Dirac operator spectrum. But it is interesting to zoom in on the microscopic behavior 

2 Many other details do not match at all. For instance the "phase transition" in the model ( |l5| ) is continuous, in 
contrast to the finite-temperature phase transition in the massless Nf = 4 theory, which is believed to be of first order. 

3 Similar numerical studies have been performed by K. Splittorff, M.Sc. thesis, The Niels Bohr Institute 1999 (unpub- 
lished) . Simulations of the macroscopic spectral density in this model can also be found in the original paper by Jackson 
and Verbaarschot |14|. 
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Figure 12: A blowup of the graphs in Fig. 11. The microscopic spectrum changes smoothly from 
Bessel-type behavior to Airy-type behavior, going through an intermediate "critical" distribution at 
T = 1 corresponding to a macroscopic power-law behavior of type p(X) ~ A 1//3 . The exact analytical 
curve for the microscopic spectral density at T = 1 has been computed in Ref. IT! 



of the Random Matrix spectral density at the soft edge of the gap. In this way, by blowing up the 



scale of the smallest eigenvalues, we obtain the plots in Fig. 12 for the same parameter values of T as 
above. One sees clearly how the universal Bessel-kernel behavior below T c turns into the also universal 
Airy- kernel above T c . We note that it is has been shown in Ref. |l4| that the massless microscopic 
spectral density of the above Random Matrix model has precisely the usual zero-temperature form, 



Ps(C(T)) 



ceo 



J Nf (((T)) 2 - J Nf -x(C(T))J Nf+1 (aT)) 



(18) 



where C(0 is simply the eigenvalues rescaled by the (T-dependent) infinite- volume spectral density 
at the origin p(0,T): 

C(T) = X2ttN P (0,T) . (19) 



In this model p(0) approaches zero with a mean-field type of behavior [14]: 



p(0,T) = p (0,0)a/1-T2 . 



(20) 



For T bigger than T c , but still close to it, we find, as expected, a deformation of the Airy-kernel. 
In fact, the microscopic behavior there smoothly interpolates between the Bessel-form and the Airy- 
form. The peaks corresponding to individual eigenvalues from the Bessel-function behavior below 
T c gradually smoothen out to become the inflection points in the spectral density of the Airy-kind 
as the soft edge moves away from the origin. To illustrate how accurately one reproduces the Airy- 
behavior in this kind of simulations, we show in Fig. [13| the soft edge prediction appropriately scaled 
to fit simulation data at T=3. The Airy-behavior is perfectly reproduced close to the edge, but with 
deviations after the first few eigenvalues (wiggles). The deviation is presumably caused by the limited 
number of eigenvalues at iV = 300, and when N is increased we expect the fit to improve. 
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Figure 13: A blowup of a T = 3 spectrum compared with the Airy-prediction. There is nice agreement 
near the edge. 

5 Conclusions 

The purpose of this study has been to find the extent to which Random Matrix Theory may be able 
to describe low-lying eigenvalue distributions and correlations between low-lying eigenvalues of the 
staggered Dirac operator at finite temperature. We have also been interested in testing the quenched 
theory with staggered fermions in the light of recent results with overlap fermions p2j , which indicated 
that the chiral finite-temperature phase transition in the quenched theory may be more subtle than 
previously expected. In the quenched case we do not see the accumulation of small Dirac operator 
eigenvalues around A = that was observed with overlap fermions. This is not entirely surprising 
in view of the insensitivity of staggered fermions to gauge field topology at these lattice couplings 
and lattice volumes ^] . With staggered fermions we do observe a strong depletion of eigenvalues 
near the origin once the temperature T reaches the pure gauge theory deconfinement phase transition 
temperature T c . This is in agreement with the conventional picture that chiral symmetry is restored 
in the quenched theory with staggered fermions at precisely the deconfinement phase transition. 

On the other hand, we find a clear difference in the behavior of the smallest Dirac operator eigenvalues 
in the quenched theory and the theory with genuine, dynamical, fermions. In the quenched case a 
small, volume independent tail of the eigenvalue distribution extends to A = while in the full theory 
the tail does not reach the origin at the couplings we have investigated. While the bulk of the 
eigenvalue distribution near the (soft) edge is roughly compatible with a square root behavior, the 
tail of small eigenvalues is considerably larger than the Airy function behavior that Random Matrix 
Theory would predict. Physically, a genuine gap in physical units in the eigenvalue spectrum above T c 
would imply the restoration of axial U(l) symmetry at these high temperatures. A very likely scenario 
is therefore that the apparent gap found with staggered fermions at finite bare coupling (3 shrinks to 
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zero in physical units as the continuum limit is reached. However, an investigation of whether this is 
indeed the case is much beyond the scope of the present paper. 
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